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We present calculations of the electronic structure of one-dimensional infinite chains and three- 
dimensional condensed matter in strong magnetic fields ranging from B = 10 12 G to 2 x 10 15 G, 
appropriate for observed magnetic neutron stars. At these field strengths, the magnetic forces on 
the electrons dominate over the Coulomb forces, and to a good approximation the electrons are 
confined to the ground Landau level. Our calculations are based on the density functional theory, 
and use a local magnetic exchange-correlation function appropriate in the strong field regime. The 
band structures of electrons in different Landau orbitals are computed self-consistently. Numerical 
results of the ground-state energies and electron work functions are given for one-dimensional chains 
Hoo, Heoo, Coc, and Feoo. Fitting formulae for the B-dependence of the energies are also provided. 
For all the field strengths considered in this paper, hydrogen, helium, and carbon chains are found 
to be bound relative to individual atoms (although for B less than a few x 10 12 G, carbon infinite 
chains are very weakly bound relative to individual atoms). Iron chains are significantly bound for 
B > 10 14 G and are weakly bound if at all at B < 10 13 G. We also study the cohesive property 
of three-dimensional condensed matter of H, He, C, and Fe at zero pressure, constructed from 
interacting chains in a body-centered tetragonal lattice. Such three-dimensional condensed matter 
is found to be bound relative to individual atoms, with the cohesive energy increasing rapidly with 
increasing B. 

PACS numbers: 31.15.Ew, 95.30.Ky, 97.10.Ld 



I. INTRODUCTION 



Youngneutron stars (ages < 10 7 years) are observed to have surface magnetic fields in the range of 10 11 -10 15 G 
[H> S S 13 > far beyond the reach of terrestrial laboratories Q . It is well known that the properties of matter can be 
drastically modified by such strong magnetic fields. The natural atomic unit for the magnetic field strength, Bq, is 
set by equating the electron cyclotron energy fiUBe = h(eB/m e c) = 11.577 Bn keV, where Bn = B/(10 12 G), to the 
characteristic atomic energy e 2 /ao = 2 x 13.6 eV (where eto is the Bohr radius): 
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B = 1^4- = 2 -3505 x 10 9 G. (1) 

For b = B I Bq > 1, the usual perturbative treatment of the magnetic effects on matter (e.g., Zeeman splitting of atomic 
energy levels) does not apply. Instead, in the transverse direction (perpendicular to the field) the Coulomb forces act 
as a perturbation to the magnetic forces, and the electrons in an atom settle into the ground Landau level. Because of 
the extreme confinement of the electrons in the transverse direction, the Coulomb force becomes much more effective 
in binding the electrons along the magnetic field direction. The atom attains a cylindrical structure. Moreover, it is 
possible for these elongated atoms to form molecular chains by covalent bonding along the field direction. Interactions 
between the linear chains can then lead to the formation of three-dimensional condensed matter 0, 0, Q ■ 

This paper is the second in a series where we present calculations of matter in strong magnetic fields using density 
functional theory. In Medin and Lai 9] (hereafter paper I), we studied various atoms and molecules in magnetic 
fields ranging from 10 12 G to 2 x 10 15 G for H, He, C, and Fe, representative of the most likely neutron star surface 
compositions. Numerical results and fitting formulae of the ground-state energies were given for Hjy (up to N — 10), 
He at (up to N = 8), Cat (up to N = 5), and Fe^ (up to N = 3), as well as for various ionized atoms. It was found 
that as B increases, molecules become increasingly more bound relative to individual atoms, and that the binding 
energy per atom in a molecule, \E^\/N, generally increases and approaches a constant value with increasing N. In 
this paper, we present density-functional-theory calculations of infinite chains of H, He, C, and Fe. Our goal is to 
obtain the cohesive energy of such one-dimensional (ID) condensed matter relative to individual atoms for a wide 
range of field strengths. We also carry out approximate calculations of the relative binding energy between ID chains 
and three-dimensional (3D) condensed matter at zero pressure. 

The cohesive property of matter in strong magnetic fields is a fundamental quantity characterizing magnetized 
neutron star surface layers, which play a key role in many neutron star processes and observed phenomena. The 
cohesive energy refers to the energy required to pull an atom out of the bulk condensed matter at zero pressure. 
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Theoretical models of pulsar and magnet ar magnetospheres depend on the cohesive properties of the surface matter 
in strong magnetic fields [hJ EH EH, El EH EH] • For example, depending on the cohesive energy of the surface matter, 
an acceleration zone ("polar gap") above the polar cap of a pulsar may or may not form, and this will affect pulsar 
radio emission and other high-energy emission processes. Also, while a hot or warm neutron star most certainly has a 
gaseous atmosphere that mediates its thermal emission, condensation of the stellar surface may occur at sufficiently 
low temperatures [1, EH]- For example, radiation from a bare condensed surface (with no atmosphere above it) has 
been invoked to exp lain the nearly perfect blackbody emission spectra observed in some nearby isolated neutron stars 
[13, EH El; HH l2l| . However, whether surface condensation actually occurs depends on the cohesive energy of the 
surface matter. 

There have been few quantitative studies of infinite chains and zero-pressure condensed matter in strong magnetic 
fields. Earlier variational calculations [22l . [23| as well as calculations based on Thomas-Fermi type statistical models 
[U HH, HH, [2?| , while useful in establishing scaling relations and providing approximate energies of the atoms and the 
condensed matter, are not adequate for obtaining reliable energy differences (cohesive energies). Quantitative results 
for the energies of infinite chains of hydrogen molecules Hoo in a wide range of field strengths (B ^> Bg) were presented 
in both Ref. [28| (using the Hartree-Fock method with the plane- wave approximation; see also Ref. || for some results 
of Heoo) and Ref. [29( (using density functional theory). For heavier elements such as C and Fe, the cohesive energies 
of ID chains have only been calculated at a few magnetic field strengths in the range of B = 10 12 -10 13 G, using 
Hartree-Fock models [3(| and density functional theory [3l[. There were discrepancies between the results of these 
works, and some (e.g., Ref. [3(|) adopted a crude treatment for the band structure (see Sec. 11110]) . An approximate 
calculation of 3D condensed matter based on density functional theory was presented in Ref. [33] • 

Our calculations of atoms and small molecule s (pape r I ) and of infinite chains and condensed matter (this paper) 
are based on density functional theory [3^, [34], l35l . [3a, H3] . In the strong field regime where the electron spins are 
aligned with each other, the Hartree-Fock method is expected to be highly accurate [3(J, HH]. However, in dealing 
with systems with many electrons, it becomes increasingly impractical as the magnetic field increases, since more 
and more Landau orbitals (even though electrons remain in the ground Landau level) are occupied and keeping track 
of the direct and exchange interactions between electrons in various orbitals becomes computational rather tedious. 
Our density-functional calculations allow us to obtain the energies of atoms and small molecules and the energy 
of condensed matter using the same method, thus providing reliable cohesive energy values for condensed surfaces 
of magnetic neutron stars, a main goal of our study. Compared to previous density-functional theory calculations 
[II [Mil 111, we use an improved exchange-correlation function appropriate for highly magnetized electron gases, 
we calibrate our density-functional code with previous results (when available) based on other methods, and (for 
calculations of condensed matter) adopt a more accurate treatment of the band structure. Moreover, our calculations 
extend to the magnetar-like field regime (B ~ 10 15 G). 

This paper is organized as follows. After briefly summarizing the approximate scaling relations for linear chains 
and condensed matter in strong magnetic fields in Sec. II, we describe our method and the basic equations in Sec. III. 
Numerical results (tables and fitting formulae) for linear chains are presented in Sec. IV. In Sec. V we describe our 
approximate calculation and results for the relative energy between ID chain and 3D condensed matter. We conclude 
in Sec. VI. Some technical details are given in the appendix. 



II. BASIC SCALING RELATIONS FOR LINEAR CHAINS AND 3D CONDENSED MATTER IN 

STRONG MAGNETIC FIELDS 



The simplest model for the linear chain is to treat it as a uniform cylinder of electrons, with ions aligned along the 
magnetic field axis. The radius of the cylinder is R and the length of a unit cell is a (which is also the atomic spacing 
along the z axis). The electrons lie in the ground Landau level, but can occupy different Landau orbitals with the 
radius of guiding center p m — (2m + l) 1 / 2 po ) where m — 0, 1, 2, ... , m max and po = (Hc/eB) 1 / 2 = (in atomic 

units). 1 The maximum Landau orbital number m max is set by Pm max = R, giving m max ~ irR 2 eB/(hc) = R 2 b/2 
(this is the Landau degeneracy in area ttR 2 ). For a uniform electron density n = Z j (itR 2 a\ the Fermi wave number 
(along z) kp is determined from n = bkpj (2ir 2 ), and the kinetic energy of the electrons in a cell is Ek = (Z/3)e' P , 
with e' F = kp/2 the Fermi kinetic energy. The total energy per atom (unit cell) in the chain can be written as [til |40| 



2ir 2 Z 3 Z^ 
a 



3b 2 R 4 a 2 



In- 



2a 
~R 



(2) 



Unless otherwise specified, we use atomic units, in which the length in ao (Bohr radius), mass in m e , energy in e 2 /ao = 2 Ry, and 
magnetic field strength in units of Bq. 
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where 7 = 0.5772 ... is Euler's constant. In Eq. ([2]), the first term is the electron kinetic energy £7. and the second 
term is the (direct) Coulomb energy (the Madelung energy for the one-dimensional uniform lattice). Minimizing 
with respect to R and a gives 

R = 1.65 Z 1/S b~ 2/5 , a/R= 2.14, 

£oo = -0.354 Z 9/5 6 2/5 . (3) 

Note that the energy @ can be written as -Eoo = —ZVq + (Z/3)e' F , where Vq is the depth of the potential well relative 
to the continuum. In equilibrium Eoo = —5Ek = — (5/3)Ze' F , and thus Vo — 2e' F . The Fermi level energy of the 
electrons in the chain relative to the continuum is then e F = e' F — Vo = —£ F = 3-Eoo/(5Z), i.e., 

e F (lD) = -0.212 Z 4/5 b 2/5 a.u. = -65.1 Z 4 ' 5 B% 5 eV. (4) 

Alternatively, if we identify the number of electrons in a cell, N e , as an independent variable, we find R = 
1.65(jV2/Z)V5&-2/5 and E x = -0.354 (Z 2 A e ) 3 / 5 6 2 / 5 . The chemical potential (which includes potential energy) 
of electrons in the chain is simply /i = ep = dE x j 'dN e , in agreement with Eq. ([!]). The electron work function is 
W=\e F \. 

A linear ID chain naturally attracts neighboring chains through the quadrupole-quadrupole interaction. By placing 
parallel chains close together (with spacing of order 6~ 2 / 5 ), we obtain three-dimensional condensed matter (e.g., a 
body-centered tetragonal lattice) |40j . 

The binding energy of the 3D condensed matter at zero pressure can be estimated using the uniform electron 
gas model. Consider a Wigner-Seitz cell with radius r$ = Z 1 ' 3 r s (r s is the mean electron spacing); the mean 
number density of electrons is n = Zj (47rr 3 /3). When the Fermi energy p F /(2m e ) is less than the electron cyclotron 
energy huiBe, or when the electron number density satisfies n < ns — {\/2i: 2 Pq y 1 = 0.0716 6 3 / 2 ( or r, > J'iB = 
1.49 Z 1 / 3 ^ -1 / 2 ), the electrons only occupy the ground Landau level. The energy per cell can be written 

3vr 2 Z 3 0.9Z 2 

^ r <) = 8653- - — ■ ^ 

where the first term is the kinetic energy and the second term is the Coulomb energy. For a zero-pressure condensed 
matter, we require dE s /dri = 0, and the equilibrium ri and energy are then given by 

n -i.qoz 1 ^- 2 / 5 , (6) 

E s ~ -0.395 Z 9/5 b 2/5 . (7) 
The corresponding zero-pressure condensation density is 

p s ~ 561 A Z- 3/5 B% 5 g cm" 3 . (8) 

The electron Fermi level energy is 

e F (3D) = ^-E s = -0.237 Z 4/5 6 2/5 a.u. = -72.7 Z^ 5 B% 5 eV. (9) 

The uniform electron gas model can be improved by incorporating the Coulomb exchange energy and Thomas-Fermi 
correction due to nonuniformity of the electron gas [1, |4l| . 

Although the simple uniform electron gas model and its Thomas-Fermi type extensions give a reasonable estimate 
for the binding energy for the condensed state, they are not adequate for determining the cohesive property of the 
condensed matter. Also, as we shall see, Eq. (|4]) or Eq. ([9]) does not give a good scaling relation for the electron 
work function when detailed electron energy levels (bands) in the condensed matter are taken into account. The 
cohesive energy Q s = E a — E s is the difference between the atomic ground-state energy E a and the condensed matter 
energy per cell E s . In principle, a three-dimensional electronic band structure calculation is needed to solve this 
problem. However, for sufficiently strong magnetic fields, such that ao/Z ^> \/2Z + lpo or B12 ^> 100 (Z/26) 3 , a 
linear ID chain is expected to be strongly bound relative to individual atoms (i.e., the cohesive energy of the chain, 
Qoc — E a — Eqo, is significantly positive) [8j. For such strong fields, the binding of 3D condensed matter results 
mainly from the covalcnt bond along the magnetic axis, rather than from chain-chain interactions; in another word, 
the energy difference |A£J S | — \E S — Eoo\ is small compared to (Joo- In the magnetic field regime where Qao is small 
or even negative, chain-chain interactions are important in deciding whether 3D condensed matter is bound relative 
to individual atoms. In this paper we will concentrate on calculating Eoo and Qoo for linear chains (Sec. Ill and 
Sec. IV). In Sec. V we shall quantify the magnitude of AE S for different elements and field strengths. 
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III. DENSITY-FUNCTIONAL-THEORY CALCULATIONS OF ID CHAINS: METHODS AND 

EQUATIONS 

Our calculations of ID infinite chains are based on density functional theory, which is well established in the strong 
magnetic field regime (B 3> -Bo) of interest here [35l |36|. Extensive comparisons of our density-functional-theory 
results for atoms and finite molecules with previous results (when available) based on different methods were given 
in paper I Q; such comparisons established the validity and calibrate the systematic error of our approach. As we 
discuss below, for infinite chains considered in the present paper, it is important to calculate the band structure of 
electrons (for different Landau orbitals) self-consistently, rather than using certain approximate ansatze as adopted 
in some previous works [301 ] . 



A. Basic equations and concepts 



Our calculations will be based on the "adiabatic approximation," in which all electrons are assumed to lie in the 
ground Landau level. For elements with nuclear charge number Z, this is an excellent approximation for 6 > Z 2 . 
Even under the more relaxed condition b 3> Z 4 / 3 , this approximation is expected to yield a reasonable total energy 
and accurate results for the energy difference between different electronic systems (atoms and chains) (see paper I). 
Also, we use nonrelativisitc quantum mechanics in our calculations, even when TiuiBe ^ m e c 2 or B > Bq = B /a 2 — 
4.414 x 10 13 G. As discussed in paper I, this is accurate as long as the electrons stay in the ground Landau level. 

In a ID chain, the ions form a periodic lattice along the magnetic field axis. The number of cells ("atoms") in 
the chain is N — > oo and the ions are equally spaced with lattice spacing a. In the adiabatic approximation, the 
one-electron wave function ("orbital") can be separated into a transverse (perpendicular to the external magnetic 
field) component and a longitudinal (along the magnetic field) component: 

*W(r) = -^=W m (r ± )f ml/k (z) . (10) 
v N 



Here W m is the ground-state Landau wave function [42| given by 



p V2^.\V2po) XP Upo 



W m (r ± ) = 7 === ( ) exp ( -JL. ) expHm^) , (11) 



which is normalized as /d 2 rj_|W m | 2 = 1. The longitudinal wave function f mvk must be solved numerically, and we 
choose to normalize it over a unit cell of the lattice: 

a/2 

\f m „ k (z)\ 2 dz = 1, (12) 

-a/2 

so that normalization of 'i'mvk is Jd 3 r \ty m vk\ 2 — 1 (here and henceforth, the general integral sign Jd 3 r refers to 
integration over the whole chain, with z from —Na/2 to Na/2). The index v = 0, 1, 2, . . . labels the different bands of 
the electron (see below) , rather than the number of nodes in the longitudinal wave function as in the atom or molecule 
case. 

The quantum number k is not present for atoms or finite molecules, but enters here because of the periodic nature 
of the electrons in the longitudinal direction. By Bloch's theorem, the electrons satisfy the periodicity condition 

f m uk(z + a) = e tka f muk (z), (13) 

and k is the Bloch wave number. Note that the longitudinal wave functions are periodic in k with period Afc = 27r/a; 
i.e., fmv,k+K{z) — fmvk{z) with K being any reciprocal vector (number, in one dimension) of the lattice, K = 2im/a 
(n is an integer). Because of this, to ensure that each wave function fmvk is unique, we restrict k to the first Brillouin 
zone, k G [—it/a, n/a]. The electrons fill each (mv) band, with spacing Afc = n/(Na), and thus the maximum number 
of electrons in a given band is N (out of the total ZN electrons in the chain). In another word, the number of 
electrons per unit cell in each (mv) band is a mu < 1 (see Sec. IIIIB[) . 
The density distribution of electrons in the chain is given by 

n(r) = \*m„k(r)\ 2 = ^ E \ W ™\ 2 (p) [ dk \f mlJk (z)\ 2 , (14) 

mvk mv 
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where the sum/integral is over all electron states, each electron occupying an (mvk) orbital. The notation | W m \ 2 (p) — 
|Wm( r J_)| 2 is used here because W m is a function of p and </> but |W m | 2 is a function of p only. The notation Jj in 
the k integral refers to the fact that the region of integration depends on the (mv) level; we will discuss this interval 
and electron occupations in Sec. IIIIBI To simplify the appearance of the electron density expression, we define the 
function 



f m Az) = ^ j i dk\f mvk {z)\\ (15) 

so that 

n{v) = Y J \W m \ 2 {p)JL{z)- (16) 

m v 

In an external magnetic field, the Hamiltonian of a free electron is 



— P + -A +- a z , (17) 

Tl e V c / 2m e c 



H 

2m 

where A = x r is the vector potential of the external magnetic field and a z is the z-component Pauli spin matrix. 
For electrons in Landau levels, with their spins aligned parallel/antiparallel to the magnetic field, the Hamiltonian 
becomes 

H=^+(n L + ^n. Be ± 1 -fu, Be , (18) 

where ul — 0, 1,2, . . . is the Landau level index; for electrons in the ground Landau level, with their spins aligned 
antiparallel to the magnetic field (so nj, = and a z — > —1), 

P 2 

H = t-^— • (19) 
2m e V 1 

The total Hamiltonian for the atom or molecule then becomes 

ff = E|^ + ^ ( 20 ) 

i 

where the sum is over all electrons and V is the total potential energy of the atom or molecule. 

In the density functional formalism, the total energy per cell of the chain is expressed as a functional of the total 
electron density n(r): 

E[n] = E K [n] + E eZ [n] + E dil [n] + E cxc [n] + E zz [n] . (21) 

Here Ek [n] is the kinetic energy of the system of non- interacting electrons, and E eZ , E<n r and E zz are the electron-ion 
Coulomb energy, the direct electron-electron interaction energy and the ion-ion interaction energy, respectively: 

N/2 

E eZ [n} = -Y Ze 2 dv^^, (22) 



j=-N/2 



\z\<a/2 



E, lu \n] = -ll drdr' ^^l , (23) 



2 



\z\<a/2 



E zz [n]=Y / — - (24) 



ja 

3=1 J 



The location of the ions in the above equations is represented by the set {zj}, with 

z j= jaz, j = (-N/2),(-N/2+l),...,0,...,N/2. (25) 
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The term E exc represents the exchange-correlation energy. In the local approximation, 

E exc [n} = / drn(r)e cxc (n) , 

J\z\<a/2 



(26) 



where £ exc (n) = e cx (n) + e corr (n) is the exchange and correlation energy per electron in a uniform electr on g as of 
density n. For electrons in the ground Landau level, the (Hartree-Fock) exchange energy can be written as [43| 



e ex («) = -Tre p Q nF(t) , 



where the dimensionless function F(t) is 

F{t) = 4 / dx 
Jo 



tan 



1 # 1 \ x I 1 



and 



n 

n B 



2n 4 p 6 n 2 , 



(27) 



(28) 



(29) 



[tib = (v / 2t'" 2 Po) 1 is the density above which the higher Landau levels start to be filled in a uniform electron gas] 
For small t, F(t) can be expanded as (4lj 



It ( 13 \ /gy 

F(i) ~3- 7 -ln4t+- f — -7-ln4tj +— ( — - 7 -ln4t) +0(t 3 \nt), 



(30) 



where 7 = 0.5772 ... is Euler's constant. We have found that the condition t <C 1 is well satisfied everywhere for 
almost all infinite chains in our calculations. The notable exceptions are the carbon chains at B = 10 12 G and the 
iron chains at B < 10 13 G, which have t ^ 1 near the center of each cell. These chains are expected to have higher t 
values than the other chains in our calculations, as they have large Z and low B 2 . 

The correlation energy of uniform electron gas in strong magnetic fields has not be calculated in general, except in 
the regime t<l and Fermi wave number kp — 2n 2 p\n 3> 1 [orn> (27r 3 pgao) _1 ]. Skudlarski and Vignale |44| use 
the random-phase approximation to find a numerical fit for the correlation energy in this regime (see also Ref. [451]): 



£corr = -_ [0.595(V&) 1/8 (1 - 1.009i 1/8 )] 
Po 



(31) 



In the absence of an "exact" correlation energy density we employ this strong-ficld-limit expression. Fortunately, 
because we are concerned mostly with finding the energy difference between atoms and chains, the correlation energy 
term does not have to be exact. The presence or the form of the correlation term has a modest effect on the atomic 
and chain energies calculated but has very little effect on the energy difference between them (see paper I for more 
details on various forms of the correlation energy and comparisons). 

Variation of the total energy with respect to the electron density, 6E[n]/6n = 0, leads to the Kohn-Sham equation: 



2m, 



■V 2 + V eS (r) 



where 



^eff(r) 



^ 2 Ze 2 



j=-N/2 1 3 



+ e 2 [ dr' t~~~jt + Mexc(n), 
r-z,- ./ r-r' 



with 



Mexc(n) 



d(ne cx 
dn 



(32) 



(33) 



(34) 



2 For the uniform gas model, t oc Z 4 / 5 _B 3 / 5 . 
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Averaging the Kohn-Sham equation over the transverse wave function yields a set of one-dimensional equations: 

h 2 d 2 



2m P dz 2 



+ V eS {z) 



fmvk(z) = £mv(k)fmvk(z) • (35) 



where 



' ~ ' \W m \ 2 { P ) , , //\ ,,\W m \ 2 (p)n(r>) 



V cS {z) = -Ze 2 £ Jdr ± \™ m }^ +e 2 jj dv^dv' 



j=-N/2 



r - r' 



+ J dr ± \W m \ 2 (p)^ c (n). (36) 

This set of equations are solved self-consistently to find the eigenvalue e mv (k) and the longitudinal wave function 
fmiik(z) for each orbital occupied by the electrons. Once these are known, the total energy per cell of the infinite 
chain can be calculated using 



Eoo = ±Y,j I dke m Ak)~^IJ z<a2 drdr' niV,nlr ' 1 



\<a/2 

f ^ Z 2 e 2 

+ / dr n(r)[e cxc (n) - p cxc (n)} + — — , (37) 

J\z\<a/2 ~[ J a 

where the interval I mu is the same as in the electron density expression, Eq. (|16p . 

Note that the electron-ion, direct electron-electron, and ion-ion interaction energy terms given above formally 
diverge for N — > oo. These terms must be properly combined to yield a finite net potential energy. Note that for an 
electron in the "primary" unit cell (—a/2 < z < a/2), the potential generated by a distant cell (centered at Zj = ja) 
can be well approximated by the quadrupole potential: 

V Q (p,z;ja) = ^^(2z 2 -p 2 ), (38) 

where Q zz is the quadrupole moment of a unit cell 

dr (2z 2 - p 2 ) n(p, z) . (39) 



|<a/2 

The Coulomb (quadrupole-quadrupole) energy between the primary cell and the distant cell is simply 



E W (>): / drn(r)V Q (p,z;ja) = ^-^. (40) 

\z\<a/2 L \3 a \ 



In our calculations, we treat distant cells with \j\ > Nq using the quadrupole approximation, while treating the 
nearby cells (|j| < Nq) exactly. Thus the (averaged) effective potential, Eq. (|36|) . becomes 



X ^ 'W m \ 2 { P ) , , ff ^ ^,\W m \ 2 {p)n(v>) 



V eS (z) = -Ze 2 £ I *x + * I i ^ dr ' 

jj^rj \ r ~ z i\ J J\z>\<a(N Q +i/2) |r-r 



+ [dr ± \W m \ 2 (p) Pexc (n)+ ( £ 4)— ?* 

\j=N Q + l 3 J a 



dr ± \W m \ 2 (p)(2z 2 - P 2 ). (41) 



The total energy per unit cell [see Eq. (37])] is given by 

e 2 f f , , , n(v)n{r') 



Eoo = 77- Y] / dke mv {k)~^— [ [ dvdv 

Zn t^Jl m „ 1 J J\z\<a/2,\z'\<a{N Q +\/2) 



r - r' 



, N Q Z 2 £ 2 ( 

/ drn(r)[e 0XC (n)-/i cxc (n)]+V— 

J\z\<a/2 y._ 



t 4 m 
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In practice, we have found that accurate results are obtained for the energy of the chain even with Nq — 1 (i.e., only 
the primary cell and its nearest neighbors are treated exactly and more distant cells are treated using quadrupole 
approximation) . 

Details of our method used in computing the various integrals above and solving the Kohn-Sham equations self- 
consistently are given in the Appendix. 



B. The electron band structure shape and occupations 



As discussed above, the electron orbitals in the chain are specified by three quantum numbers: m, v, k. While m, v 
are discrete, k is continuous. In the ground state, the electrons will occupy the (mvk) orbitals with the lowest energy 
eigenvalues e mv (k). To determine the electron occupations and the total chain energy, it is necessary to calculate 
the e m v{k) energy curves. Here we discuss the qualitative property of these energy curves (i.e., the electron band 
structure) using the theory of one-dimensional periodic potentials (see, e.g., Ref. p6|). 

Like the wave functions, the energy curves are periodic, with £ m „(fc + K) = e mv (k), where K is lix ja multiplied by 
any integer. The energy curves are also symmetric about the Bragg "planes" ("points" in ID) of the reciprocal lattice, 
e mv (K — k) = e rnu (k). Thus we can determine the entire band structure of the electrons by calculating it between 
any two Bragg points. Since we have chosen to limit our calculation to the first Brillouin zone k G [—ir/a,w/a], we 
only need to consider the domain k G [0, 7r/a]. 

For a given m, the energy curves lie in bands which do not overlap and increase in energy with increasing v (see 
Fig. [I]). These bands are bounded by the energy values at the Bragg points, such that in each band the energy 
increases/decreases monotonically between the two points. The direction of this growth alternates with v. For the 
v = band, the energy is at a minimum for k = and increases to a local maximum at k — 7r/a; for the v = 1 band, 
the energy curve is at a minimum for k = it / a and grows to a maximum at k = 0, etc. These properties are depicted 
in Fig. [0 

Also shown in the figure is the Fermi level energy ep of the electrons in the infinite chain. The electrons occupy 
all orbitals {mvk) with energy less than ep. For each (mv) band, we define the occupation parameter <x m „, which 
gives the number of electrons that occupy this band per unit cell [i.e., the number of electrons that occupy the (mv) 
band in the whole chain is a mv N]. Since the maximum possible number of electrons in each (mv) band is N, we have 
(r w < 1. Because there are ZN electrons total in the chain, these occupation numbers are subject to the constraint 



(43) 



It is also useful to define for each (mv) level the Fermi wave number h™", such that the electrons fill up all allowed 
orbitals between the minimum-energy Bragg point (k — for even v and k = ix j a for odd v) and hp" . The occupied 
fc's are therefore 
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(44) 



(45) 



for odd v, plus the corresponding reflection about the Bragg point k = 0. For a completely filled band (as illustrated 
in Fig. [T]for the v = band), a mv — 1 and kp 1 " = ir/a (for v = even) or (for v — odd); for a partially filled band 
(the v = 1 band in Fig. [1]), 

s m v(k r r) = £f • (46) 
With the allowed k values specified, the k integration domain in Eqs. (fT4"|) , (|15[) . ([57)) and (|4"2|) is given by 



dk 



2 f Q F dk, v even, 



odd. 



(47) 



Note that the Fermi level energy ep and various occupation numbers er TO „ must be calculated self-consistently. In 
principle, they should be determined by minimizing the total energy with respect to <j mv subject to the constraint 
Eq. g3|), i.e., 



E[n; a mv ] - e F 



Z 



0. 



(48) 
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FIG. 1: A schematic diagram showing the electron band structure for a particular m value. In this example, the first band 
(y = 0) is fully occupied (cr TO o = 1) while the second band (y — 1) is partially filled {<r m i < !)• 



Since 



dn(r) 7T dnjr) 2 2 



(49) 



Eq. ESJ yields 



2m P g?z 2 



(50) 



Comparing this to Eq. ([55]) . we find e mv (k™ v ) = Ef, which is Eq. PS)) . This shows that using Eq. (|4^|) to find £f 
minimizes the total energy of the system. 
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Az=a 



FIG. 2: A schematic diagram showing the shapes of the longitudinal wave functions of electrons in different bands at k = 
and k = n/a. 



C. The complex longitudinal wave functions 

The longitudinal electron wave function f mv k(z) satisfies the Kohn-Sham equations (|35p subject to the periodicity 
condition Eq. (| 13[> . or equivalently, the cell boundary condition 

f m uk(a/2) = e* ka f mvk (-a/2). (51) 

Since the electron density distribution n(r) is periodic across each cell and symmetric about each ion, the following 
boundary condition is also useful: 

\fmvk(z)\'\z=a/2 = \fm V k{z)\'\z=-a/2 = 0. (52) 

Due to the complex boundary condition Eq. (|5ip . the wave function f m vk is complex for general fc's. The exceptions 
are k = and fc = n/a: For k = 0, the boundary condition becomes f m vk(a/2) = fmvk(—a/2), and for k = ir/a we 
have /mi/fc(«/2) = — ,fmvk{— a/2). Thus for k — and 7r/a, we can choose the longitudinal wave functions to be real. 
The general shapes of these wave functions (for different bands) are sketched out in Fig. [2] We see that at the Bragg 
points, between the two states with the same number of nodes, the one that is more concentrated near the ion has 
lower energy than the other state; this difference gives rise to the band gap. The k — 0,ir/a eigenvalues e mu and 
eigenfunctions can be calculated in the domain < z < a/2 with the boundary condition / mj/ fc(0) = or f' mvk {$) = 0. 

The electron wave functions for general fc's are more difficult to compute as they have complex boundary conditions. 
Our procedure for calculating these wave functions and their corresponding electron energies is as follows: For each 
energy band (mi'), the electron eigenstates at k — and k — n/a are first found (see above). For every energy between 
£mv(k = 0) and e mv (k = -k /a), we find the wave function that solves the Kohn-Sham equation while satisfying the 
symmetric/periodic density condition Eq. (I52p . More precisely, we choose / = 1 (up to a normalization constant) and 
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guess f — ig (where g is a real number) at z = a/2 (thus |/|' = is satisfied at z = a/2), and then integrate the 
Kohn-Sham equation to z = —a/2; we adjust g so that |/|' = is satisfied at z = —a/2. Once the wave function is 
obtained, we determine its k value from the Bloch boundary condition Eq. (|5ip . Through this method we find s mv (k) 
as a function of k for each (mv) band. 

Some examples of our computed s mv (k) are shown in Figs. [3] and [4] To show that our calculations are consistent 
with theoretical models, we have included several model fits for the electron energy curves: the tight-binding fit in 
Fig. [31 which has the form 

e mv {k) ~ ci + c 2 cos(fca) (53) 
[see Ref. |46[, Eq. (10.19)], and the weak-periodic-potential fit in Fig. 0] which has the form 

e mv (k) ~ ci + l -[k 2 /2 + (2n/a - kf /2] - l -{[{2n/a - k) 2 /2 - k 2 /2} 2 + c 2 } 1 ' 2 (54) 

[see Ref. [46} . Eq. (9.26)]. The constants ci and c 2 in the formulas are fit to the two endpoints of the energy curves, 
£mv(0) and £ mv (7r/a). The tight-binding model fits well for the most-tightly-bound electron bands in our calculations, 
while the weak-periodic-potential model fits well for all of the other bands. Note that for k <C 7r/a, the electron energy 
can be approximately fit by £ m „(fc) = e mv (0) + k 2 /2, as would be the case if the wave functions were of the form 
f m v{z)e lkz — this is the ansatz adopted by Neuhauser et al. [3(| in their Hartree-Fock calculations. But obviously 
for larger k, this is a rather bad approximation. We suggest that approximate treatment in the band structure may 
account for a large part of the discrepancies among cohesive energy results in previous works. For example, the 
disagreement between Ref. [3l[ [where e m „(fc) was calculated for a few values of k and then fit to a simple expression] 
and Ref. [301 ] (where a k 2 dependence for the electron energy was assumed) on whether or not carbon is bound at 
B\2 = 5 is due to the band structure model, not to the fact the former used the density functional theory while the 
latter used the Hartree-Fock method. 



IV. RESULTS: ONE-DIMENSIONAL CHAINS 



In this section we present our results for hydrogen, helium, carbon, and iron infinite chains at various magnetic 
field strengths between B — 10 12 G and 2 x 10 15 G. For each chain, data is given in tabular form for the ground-state 
energy (per unit cell) B^, the equilibrium ion separation a, and the electron Fermi level energy Sp (the electron work 
function is W — \s F \)- We provide relevant information for the electron occupations in different bands, such as the 
number of Landau orbitals and the number of fully occupied bands (see below for specific elements). We also give 
the ground-state energy of the corresponding atom, E a , so that the cohesive energy of each chain can be obtained, 
Qoo = B a Eoq. 

For each chain and atom we provide numerical scaling relations for the ground-state energy and Fermi level energy 
as a function of the magnetic field, in the form of scaling exponents /3 and 7, with 

E a , E x cx B 13 , e F cx EP . (55) 

We also give the rescaled, dimensionless energy Ego, and equilibrium ion separation a defined by [see Eq. ([3"|1] 

Eoo ~ E x Z g / 5 b 2/5 a.u., a ~ a Z 1 ' 5 ^ 2 ' 5 a.u.. (56) 

We shall see that the scaling relations in Eq. (|56|) with E^ ~ const, and a ~ const, represent a reasonable 
approximation to our numerical results, although such scaling formulae are not accurate enough for calculating the 
cohesive energy Qoa = Ea — E^. However, Eq. ^ or Eq. © for the Fermi level energy based on the uniform gas model 
is not a good representation of our numerical results. In paper I [9j we have shown that as N increases, the energy 
per atom in the Hjv (or He^r, Cat, Fe^y) molecule, En/N, gradually approaches a constant value. The infinite chain 
ground-state energy E^ found in the present paper is consistent with the large- N molecule ground-state energy limit 
En/N obtained in paper I (see the related figures in the following subsections). Since finite molecules and infinite 
chains involve completely different treatments of the electron states, the consistency of Eqo and E^/N provides an 
important check of the validity of our calculations. 

Other comparisons can be made between the infinite chains and finite molecules. For example, our results of ion 
separation a and scaling constant (3 are consistent between infinite chains and finite molecules. Also, we find that if 
the isolated atom has electrons in v = and v = 1 orbitals, then the corresponding infinite chain will have electrons 
in v = and v — 1 bands; if the isolated atom only has electrons in v = orbitals, the corresponding infinite chain 
will have electrons only in v = bands. 

We have compared our cohesive energy results with those of other work, whenever available. These comparisons 
are presented in the following subsections. 
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FIG. 3: The electron energy of the (m,u) = (0,0) band for the carbon infinite chain at B12 = 1. The tight-binding model fit 
for this level is shown as a dashed line [see Eq. ([53}], and the dotted line shows the free electron result eoo{k) — £oo(0) = fc 2 /2. 



A. Hydrogen 

Our numerical results for H are given in Table [H Examples of the energy curves of various Hjv molecules and Hoc 
at B12 — 1 are depicted in Fig. [5j The minimum of each energy curve determines the equilibrium ion separation 
in the molecule/chain. Figure [6] compares the molecular and infinite chain energies at various field strengths, and 
shows that as N increases, the energy per atom in the H^v molecule asymptotes to -Eoo- Figure [7] gives the occupation 
number er m o of different Landau orbitals at various field strengths. Only the v = bands are occupied, none of these 
are completely filled (cr m o < l)j an d the v >1 bands are empty (er TO i = 0). We see that as B increases, the electrons 
spread into more Landau orbitals, thus the number of m states occupied by the electrons (n m in Table [J) increases. 
Approximately, since the chain radius R oc &~ 2 / 5 and R ~ (2n m — l) 1 / 2 /^ 1 / 2 (the electrons occupy the Landau orbitals 
with m = 0, 1, 2, . . . , n m ~ 1), we have n m oc b 1 ^ 5 . Table [J shows that for B12 ^ 10 our results for and a are well 
fit by 

Eoo -529 B% 374 eV, a = 0.091 B^ A0 a (57) 

[where -B13 = i?/(10 13 G)], similar to the scaling of Eq. (|5S|) . The electron work function W = \sp\ does not scale as 
Eq. 0, but is a fraction of the ionization energy of the H atom, \E a \. Note that \E a \ is not well fit by a power law 
(oc B 13 ), but is well described by \E a \ oc (ln6) 2 (accurate fitting formulae for \E a \ are given in, e.g., Ref. [471]). 

At B 12 = 1, 10, 100, we find cohesive energies of Q x = E a - E x = 59.6, 219.7, 712.7 eV (see Table 0. At those 
same fields, Lai et al. 2g] find cohesive energies of 28.9, 141, 520 eV. At B12 = 0.94, Relovsky and Ruder [2^] 
find a cohesive energy of 47.1 eV. We expect our H calculation (and that of Ref. (29[) to overestimate the cohesive 
energy since an exchange-correlation functional is used in the chain calculation while none is required for the H atom. 
But we also expect the result obtained in Ref. [28[ to somewhat underestimate the cohesive energy since a uniform 
(longitudinal) electron density was assumed. 
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FIG. 4: The electron energy of the (m, v) = (0, 0) band for the iron infinite chain at B12 = 2000. The weak-periodic- 
potential model fit for this level is shown as a dashed line [see Eq. (|54[) ]. and the dotted line shows the free electron result 
£oo(fc) -eoo(0) = k 2 /2. 



TABLE I: The ground-state energy (per unit cell) (in units of eV), ion separation a (in units of Bohr radius 00), the 
number of occupied Landau levels n m , and the Fermi level energy ef (in eV) of ID infinite chains of hydrogen, over a range 
of magnetic field strengths. The ground-state energy of individual hydrogen atoms, E a (in units of eV), is also provided for 
reference. The dimensionless energy Boo and ion separation a are calculated using Eq. (|56jl . The scaling exponents ft and 7, 
defined by E a , Eoo oc B , and ef oc B 1 , are calculated over the three magnetic field ranges provided in the table: B12 = 1 — 10, 
10 — 100, 100 — 1000 (the exponent in the B12 = 1 row corresponds to the fit over B12 = 1 — 10, etc.). The occupation of 
different (mu) bands is designated by the number n m : the electrons occupy Landau orbitals with m = 0,1,2, ... , n m — 1, all 
in the v — band; see Fig. [7] 



B12 


H 

E a f3 


Eoc E x (3 a a n m ef 7 


1 


-161.4 0.283 


-221.0 -0.721 0.379 0.23 2.6 6 -85.0 0.28 


10 


-309.5 0.242 


-529.2 -0.688 0.374 0.091 2.6 10 -165 0.27 


100 


-540.3 0.207 


-1253.0 -0.648 0.374 0.037 2.6 16 -311 0.26 


1000 


-869.6 


-2962 -0.610 - 0.0145 2.6 26 -571 - 



B. Helium 

Our numerical results for He are given in Table [TTJ Figure [H] compares the molecular and infinite chain energies 
at various field strengths, and shows that as N increase, the energy per atom in the He at molecule approaches E x 
for the infinite chain. Figure [9] gives occupation number er TO o of different Landau orbitals at various field strengths. 
As in the case of H, only the v — bands are occupied, and the number of Landau states required (n m in Table |Hj 
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FIG. 5: The energies (per atom or cell) of various H molecules and infinite chain as a function of ion separation a at B12 = 1. 
The results of finite molecules are based on paper I [9(. The energy of the H atom is shown as a horizontal line at —161.4 eV. 



increases with increasing B, with n m oc Z 2 / 5 ^ 1 / 5 . Table ITT1 shows that for J3 12 > 10, 

Eoo ~ -1252 B^ 382 eV, a = 0.109 B^ AQ a , (58) 

similar to the scaling of Eq. (|56|) . The electron work function W = \ef\ does not scale as Eq. (|4|, but is a fraction 
of the ionization energy: Using a Hartree-Fock code (e.g., Ref. [28| ) we find that at £>i 2 = 1,10,100,1000 the He 
atomic energies are —575.5, —1178.0, —2193, —3742 eV. The He + (i.e., once-ionized He) energies at these field 
strengths are —416.2, —846.5, —1562.0, —2638 eV. Therefore, the ionization energies of He at these field strengths 
are 159.3, 331.5, 631, and 1104 eV, respectively. 

At Byi = 1, we find a cohesive energy of 58.9 eV (see Table HI)) . At the same field, Neuhauser et al. [30] (based 
on the Hartree-Fock model) find a cohesive energy of 25 eV, and Miiller [23| (based on variational methods) gives a 
cohesive energy of 50 eV. At £> 12 = 0.94, Relovsky and Ruder [2^| (based on density functional theory) find a cohesive 
energy of 56.6 eV. At B12 = 5 Jones |3li fi nds a cohesive energy of 220 eV, which is close to our value. That our 
results agree best with those of Refs. [29|, |31| is expected, as we used a similar method to find the ground-state atomic 
and chain energies. Similar to the finite He molecules (paper I), we expect our density- functional-theory calculation 
to overestimate the cohesive energy, but we also expect the result of Ref. [3(| to underestimate Qoo. 



C. Carbon 



Our numerical results for C are given in Table IIII1 Figure [TO] compares molecular and infinite chain energies at 
various field strengths, showing that as N increase, the energy per atom in the Cjv molecule approaches Eoo for the 
infinite chain. Figure ITTI gives the occupation number a m Q of different Landau orbitals at various field strengths. As 
in the case of H and He, only the v = bands are occupied, although for C at B12 = 1, the m = and m = 1 
bands (both with v — 0) are fully occupied (thus rif = 2 in Tabic IIII[) . The number of Landau states required (n m in 
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FIG. 6: The molecular energy per atom, \En\/N, for the Hjv molecule, as a function of N at several different field strengths. 
The results of finite molecules are based on paper I |9|]. As N increases, En/N asymptotes to E^. To facilitate plotting, the 
values of \E\\ (atom) at different magnetic field strengths are normalized to the value at B\ 2 = 1, 161.4 eV. This means that 
a = 1 for B12 = l,a = 161.4/309.5 for B 12 = 10, a = 161.4/540.3 for B 12 = 100, and a = 161.4/869.6 for B 12 = 1000. 



Tabic HlT| increases with increasing B, approximately with n m oc Z 2 / 5 ^ 1 / 5 . Table UTT1 shows that for Byi ^ 10, 

Eoo ~ -10 300 B? 3 387 eV, a = 0.154 B^ A3 a . (59) 

Note that these expressions are more approximate than for H and He. The electron work function W — \ef\ does 
not scale as Eq. (U|), but is a fraction of the ionization energy: from paper I 0], the ionization energies of C at 
B u = 1, 10, 100, 1000 are 174, 430, 990, and 2120 eV, respectively. 

At B12 = 10, we find a cohesive energy of 240 eV (see Table UTTT) . At B\2 — 8.5, Relovsky and Ruder [29| give 
a cohesive energy of 240 eV. At £>i 2 = 5 Jones [3l[ finds a cohesive energy of 100 eV; at the same field (using our 
scaling relations), we find a cohesive energy of 100 eV (±30 eV). Neuhauser et al. [3(|, on the other hand, find that 
carbon is not bound at Byi — 1 or 5. This is probably due to the approximate band structure ansatz adopted in 
Ref. [3(| (see Sec lIIlCI) : for fully occupied bands, the approximation that s mv (k) increases as fc 2 /2 is invalid and can 
lead to large error in the total energy of the chain. 



D. Iron 



Our numerical results for Fe are given in Table HVl The electron density profile at various field strengths is shown 
in Figs. 1121 and [TBI As the magnetic field increases the density goes up, for two reasons. First, the equilibrium ion 
separation decreases. Second, the electrons become more tightly bound to each ion, in both the p and z directions 
(the electrons move closer to each ion faster than the ions move closer to each other). It is interesting to note that the 
peak density at a given z is not necessarily along the centeral axis of the chain (p = 0), but gradually moves outward 
with increasing z. 
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FIG. 7: The occupation numbers of each m level of hydrogen infinite chains, for various magnetic field strengths. The data 
points are plotted over the curves to show the discrete nature of the m levels. Note that only the v — bands are occupied by 
the electrons. 



TABLE II: The ground-state energy (per unit cell) Eoo (in units of eV), ion separation a (in units of Bohr radius ao), the 
number of occupied Landau levels n m , and Fermi level energy ef (in eV) of ID infinite chains of helium, over a range of 
magnetic field strengths. The ground-state energy of individual He atoms, E a (in units of eV), is also provided for reference 
(this is based on the density-functional-theory calculation of Ref. Q]). The dimensionless energy Eoo and ion separation a 
are calculated using Eq. (|56[) . The scaling exponents /3 and 7, defined by E a , E x oc B , and ef <x B 1 , are calculated over 
the three magnetic field ranges provided in the table: B12 = 1 — 10, 10 — 100, 100 — 1000 (the exponent in the B12 = 1 row 
corresponds to the fit over B12 = 1 — 10, etc.). The occupation of different (mv) bands is designated by the number n m : the 
electrons occupy Landau orbitals with m = 0, 1, 2, . . . , n m — 1, all in the v = band; see Fig. U] Note that all of the He atoms 
here also have electrons only in the v = states. 



B\2 


He 

E a 


Heoo 

Eoo Eoo ft a a rim. ef 7 


1 


-603.5 0.317 


-662.4 -0.621 0.385 0.28 2.7 9 -85.0 0.29 


10 


-1252.0 0.280 


-1608.0 -0.600 0.382 0.109 2.7 14 -167 0.27 


100 


-2385 0.248 


-3874 -0.575 0.382 0.043 2.7 23 -310 0.26 


1000 


-4222 


-9329 -0.552 - 0.0175 2.7 39 -568 - 



The energy curves for Fe2, Fe3 (calculated in paper I @), and Feoo at B i2 = 500 are shown in Fig. [14] Figure [TBI 
compares the molecular and infinite chain energies at various field strengths, showing that as N increases, the energy 
per atom in the Fe^y molecule approaches -Eoo for the infinite chain. Figure [16] gives the occupation number a ml , of 
different bands at various field strengths. For B12 > 100, only the v = bands are occupied; for such field strengths, 
the Fe atom also has all its electrons in the tightly bound v — states (see Table ITV]) . At B12 — 100, the number 
of fully occupied bands is — 7 (m = 0, 1, 2, . . . , 6, all with v = 0). As B increases, the electrons spread to more 
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FIG. 8: The molecular energy per atom, \En\/N, for the Hejv molecule, as a function of N at several different field strengths. 
The results of finite molecules are based on paper 10. As TV increases, En/N asymptotes to E^. To facilitate plotting, the 
values of \E\\ (atom) at different magnetic field strengths are normalized to the value at B\ 2 = 1, 603.5 eV. This means that 
a = 1 for B12 = l,a = 603.5/1252.0 for B 12 = 10, a = 603.5/2385 for B 12 = 100, and a = 603.5/4222 for B 12 = 1000. 



Landau orbitals, and the number of occupied m-states rim increases, approximately clS 77.777 OC Z ' 5 b ' . Note that at 
the highest field strength considered, the electrons occupy m — 0, 1, 2, . . . , 156 — keeping track of all these Landau 



orbitals (rim^ — 157) is one of the more challenging aspects of our computation. Table ITVl shows that for Byi ^ 100, 



-356 5° 4 374 keV, 



a = 0.107 B7? A3 a 



(60) 



[where B14 = _B/(10 14 G)]. These scaling expressions are more approximate than for H and He. The electron work 
function W = \sf\ does not scale as Eq. (j4]), but is a fraction of the ionization energy: from paper I ||, the ionization 
energies of Fe at B 12 = 100, 500, 1000, 2000 are 1.2, 2.5, 3.4, and 5.5 keV, respectively. 

Note that at B12 = 5 and 10, the cohesive energy (Qoo = E a — E^) of the iron chain is rather small compared 
to the absolute value of the ground-state energy of the atom (|J5 a |) or chain (lE^I). For these field strengths, our 
formal numerical result for the cohesive energy is at or smaller than the standard error of our computations (0.1% 
of \E a \ or I-EqoI), so we have redone the calculations using more grid and integration points such that the atomic 
and chain energies reported here for these field strengths are accurate to at least 0.02% of \E a \ or lE^] (see the 
Appendix). Although these more- accurate cohesive energies are (barely) larger than the error in our calculations, 
there are of course systematic errors introduced by using density functional theory which must be considered. It is 
very possible that a similar, full-band-structure calculation using Hartree-Fock theory would find no binding. In any 
case, for such "low" field strengths [Bn ^ 10) the exact result of our one-dimensional calculation is not crucial, since 
in the three-dimensional condensed matter the additional cohesion resulting from chain-chain interactions dominates 
over Qoo, as we will show in Se c. fVl 

At Byi = 5, Neuhauser et al. (30t| and Jones [3lj found that iron is not bound, while we find that it is barely bound. 
At B12 — 10, Jones [32j calculated the cohesive energy for three-dimensional condensed matter, so we compare our 
results with those of Ref. [32| in Sec. |V] We have not found any quantitative calculations of cohesive energies for iron 
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FIG. 9: The occupation numbers of each m level of infinite He chains, for various magnetic field strengths. Only the v = 
bands are occupied by the electrons. Note that for B\2 — 1, the m — 8 orbital has a rather small occupation, ago ~ 0.006; if 
£f were slightly more negative, this orbital would be completely unoccupied. 

TABLE III: The ground-state energy (per unit cell) Eoa (in units of eV), ion separation a (in units of Bohr radius ao), electron 
occupation numbers (n m ,n/), and Fermi level energy sf (in eV) of ID infinite chains of carbon, over a range of magnetic field 
strengths. The ground-state energy of individual C atoms, E a (in units of eV), is also provided for reference (this is based on 
the density-functional-theory calculation of Ref. *9|). The dimensionless energy E^ and ion separation a are calculated using 
Eq. (|56[1 . The scaling exponents (5 and 7, defined by E a , E^, oc B 13 , and £f <x B 1 , are calculated over the three magnetic 
field ranges provided in the table: B12 — 1 — 10, 10 — 100, 100 — 1000 (the exponent in the B12 = 1 row corresponds to the fit 
over B12 = 1— 10, etc.). The occupation of different (mv) bands is designated by the notation (n m ,n/): the electrons occupy 
Landau orbitals with to = 0, 1,2, ... , n m — 1, all with v — 0; the number of fully occupied (a mv = 1) bands is denoted by n/; 
see Fig. 1111 Note that all of the C atoms here also have electrons only in the v = states. 



B12 


c 

E a p 


Coo 

Eoo Eoa (3 a a (n m ,nf) e F 7 


1 


-4341 0.366 


-4367 -0.567 0.373 0.49 3.9 (12,2) -92.8 0.27 


10 


-10075 0.326 


-10315 -0.533 0.385 0.154 3.1 (23,0) -173 0.25 


100 


-21360 0.287 


-25040 -0.515 0.389 0.056 2.8 (41,0) -306 0.25 


1000 


-41330 - 


-61320 -0.502 - 0.022 2.7 (69,0) -539 - 



at field strengths larger than 10 13 G. 
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FIG. 10: The molecular energy per atom, \En\/N, for the Cjv molecule, as a function of N at several different field strengths. 
The results of finite molecules are based on paper I 0. As N increases, En/N asymptotes to i?oo- To facilitate plotting, the 
values of \E\\ (atom) at different magnetic field strengths are normalized to the value at B\ 2 = 1, 4341 eV. This means that 
a = 1 for B 12 = l,a = 4341/10075 for B 12 = 10, a = 4341/21 360 for B 12 = 100, and a = 4341/41 330 for B 12 = 1000. 



TABLE IV: The ground-state energy (per unit cell) Eaa (in units of keV), ion separation o (in units of the Bohr radius do), 
electron occupation numbers (n^npn^nj 1 '), and Fermi level energy ef (in eV) of ID infinite iron chains, over a range of 
magnetic field strengths. The ground-state energy of individual Fe atoms, E a (in units of keV), is also provided for reference 
(this is based on the density- functional-theory calculation of Ref. |J). The dimensionless energy E^ and ion separation a 
are calculated using Eq. (|56p . The scaling exponents j3 and 7, defined by E a , E x oc B B , and ef oc B" 1 , are calculated over 
the three magnetic field ranges provided in the table: B\ 2 = 1 — 10, 10 — 100, 100 — 1000 (the exponent in the B\ 2 = 1 row 
corresponds to the fit over B\ 2 = 1 — 10, etc.). For atoms the electron configuration is specified by the notation (no,ni) (with 
no + ni — Z — 26), where no is the number of electrons in the v = orbitals and ri\ is the number of electrons in the v — 1 
orbitals. For infinite chains, the occupation of different (mu) bands is designated by the notation (n$ , ; nS' , where 
ri~m is the total number of occupied v = orbitals (from m = to m = n$ — 1), and rim the corresponding number for the 
v = 1 orbitals; {nP) is the number of fully occupied (a mv = 1) v = (u = 1) orbitals. Note that for Bi 2 > 100, only the 
v = states are occupied in the Fe atom, and only the v = bands are occupied in the Fe chain; see Fig. 1161 
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0.503 
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FIG. 11: The occupation numbers of each m level of infinite C chains, for various magnetic field strengths. Only the v = 
bands are occupied by the electrons. Note that for B12 = 1, the m = and m = 1 bands are completely filled. 



V. CALCULATIONS OF THREE-DIMENSIONAL CONDENSED MATTER 



For the magnetic field strengths considered in this paper (B > 10 12 G), H and He infinite chains are significantly 
bound relative to individual atoms. Additional binding energy between 3D condensed matter and ID chain is expected 
to be small [28[ (see below). Thus the cohesive energy of the 3D condensed H or He, Q s — E a — E s (where E s is 
the energy per cell in the 3D condensed matter), is close to = E a — E x , the cohesive energy of the ID H or 
H chain. For C and Fe at relatively low magnetic fields (e.g., C at B\2 < 10 and Fe at Bn ^ 100), ID chains 
are not significantly bound relative to atoms and additional cohesion due to chain-chain interactions is important in 
determining the true cohesive energy of the 3D condensed matter. Indeed, for Fe at B\i = 5, 10, our calculations of 
ID chains give such a small (see Table ITVT) that it is somewhat ambiguous as to whether the Fe condensed matter 



is truly bound relative to individual atoms. In these cases, calculations of 3D condensed matter is crucial [321 ]. 

In this section, we present an approximate calculation of the relative binding energy between 3D condensed matter 
and ID chains, AE* — E, — E^. 



A. Method 



To form 3D condensed matter we place the infinite chains in parallel bundles along the magnetic field. We consider 
a body-centered tetragonal lattice structure; i.e., the chains are uniformly spaced in over a grid in the xy plane 
(perpendicular to the magnetic axis), with every other chain in the grid shifted by half a cell (Az — a/2) in the 
z direction. The transverse separation between two nearest neighboring chains is denoted by 2R, with R to be 
determined. 

To calculate the ground-state energy of this 3D condensed matter, we assume that the electron density calculated 
for an individual ID chain is not modified by chain-chain interactions, thus we do not solve for the full electron 
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FIG. 12: The density distribution of electrons in the iron infinite chain at four different magnetic field strengths (labeled on 
the graphs). The density is shown as a function of p for five equally spaced z points from the center of a cell (z = 0) to the 
edge of that cell (z — a/2). 



density in the 3D lattice self-consistently. In reality, for each Landau orbital the transverse wave function of an 
electron in the 3D lattice is no longer given by Eq. (Ill)) (which is centered at one particular chain) , but is given by a 
superposition of many such Landau wave functions centered at different lattice sites and satisfies the periodic (Bloch) 
boundary condition. The longitudinal wave function f m uk{z) will be similarly modified. Our calculations show that 
the equilibrium separation (2R) between chains is large enough that there is little overlap in the electron densities of 
any two chains, so we believe that our approximation is reasonable. 

Using this approximation, the electron density in the 3D lattice is simply the sum of individual infinite chain 
electron densities: 

n 3 D(r) = ^n(r-r#), (61) 

ij 

where n(r) are the electron density in the ID chain (as calculated in Sees. IIIlllVI) . the sum over ij spans all positive 
and negative integers, and 

r y = 2Mx + 2Ejy+|[i,j]z (62) 

represents the location of the origin of each chain (the notation — 1 when i + j = odd, and = when i + j = 
even). In practice, the chain-chain overlap is so small that we only need to consider neighboring chains. The density 
at a point in the positive xyz octant of a 3D unit cell is approximately given by 

7i 3D (r) ~ n(r) + n(r - 2i?x - o/2z) + n(r - 2Ry - a/2z) + n(r - 2i?x - 2Ry) . (63) 

The energy (per unit cell) AE^d (R) of the 3D condensed matter relative to the ID chain consists of the chain-chain 
interaction Coulomb energy Ai?coui and the additional electron kinetic energy AEk and exchange-correlation energy 
AE CXC due to the (slight) overlap of different chains. The dominant contribution to the Coulomb energy comes from 
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FIG. 13: The density distribution of electrons in the iron infinite chain at four different magnetic field strengths (labeled on 
the graphs). The density is shown as a function of z for five equally spaced p points from the center of a cell (p — 0) to the 
guiding center radius of the highest occupied m level (p = pm max ). The p points are given in units of ciq. 



the interaction between nearest-neighboring cells. For a given cell in the matter, each of the eight nearest-neighboring 
cells contributes an interaction energy of 

E nn = E e z,nn + -Edir,nn + EzZ.nn , (64) 

where 

E eZ , nn = -Ze 2 [ dr n(r) , (65) 
J\z\<a/2 |r-r nn | 

2 J J\ z \<a/2,\z>\<a/2 |r- (r' + r nn )| 



_1ZV_1 ZV 
EzZ ' nn ~ 2 - V(a/2)H(2ffl ' (6?) 

and r nn is the location of the ion in a nearest-neighboring cell, for example 

r nn = 2i?x+^z. (68) 

More distant cells contribute to the Coulomb energy through their quadrupole moments. The classical quadrupole- 
quadrupole interaction energy between two cells separated by a distance d is 

Eqq (d,6) = ^0^(3-30 cos 2 9 + 35 cos 4 9) , (69) 
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FIG. 14: The energy per cell as a function of the ion separation for an infinite Fe chain at B12 = 500. The molecular energy 
per atom versus ion separation for the Fe2 and Fe3 molecules at the same field strength (based on calculations in paper I) are 
also shown. The energy of the Fe atom is shown as a horizontal line at —637.8 keV. 



where Q zz is given by Eq. (|39[) and 6 is the angle between the line joining the two quadrupoles and the z axis. The 
total contribution from all nonneighboring cells to the Coulomb energy is then 

\ E QQ (r ijk ), (70) 

(ijk) 

where 

Tij k =Tij + akz, d=\r ijk \, cos9 = - ^-Jj^— ~ , ( 71 ) 

and the sum in Eq. (|70[) spans over all positive and negative integers except those corresponding to the nearest 
neighbors. 

In the density functional theory, the kinetic and exchange-correlation energies depend entirely on the electron 
density. These energies differ in the 3D condensed matter from the ID chain because the overall electron density 
H3d(f) [see Eq. (|6ip ] within each 3D cell is (slightly) larger than n(r) due to the overlap of the infinite chains. Since 
we do not solve for the electron density in the 3D condensed matter self-consistently, we calculate the kinetic energy 
difference using the local (Thomas-Fermi) approximation: 

AE K (R)= drn m (r)s K (n 3V )- drn(r)e K {n). (72) 

J\z\<a/2;\x\,\y\<R. J\z\<a/2 

Here exip) is the (Thomas- Fermi) kinetic energy (per electron) for an electron gas at density n, and is given by (e.g., 
1) 

e,w ^ 2(2 ; w = f b ^t, (73) 

6m e Spo 
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FIG. 15: The molecular energy per atom, |i?jv|/iV, for the Feiv molecule, as a function of N at several different field strengths. 
The results of finite molecules are based on paper I |9|]. As N increases, En/N asymptotes to Eoo. To facilitate plotting, the 
values of (atom) at different magnetic field strengths are normalized to the value at Bi 2 = 100, 354.0 keV. This means that 
a = 1 for B 12 = 100, a = 354.0/637.8 for B 12 = 500, a = 354.0/810.6 for B 12 = 1000, and a = 354.0/1021.5 for B 12 = 2000. 



where t is given by Eq. (|29p . Note that the regions of integration in the xy direction are different for the two terms 
in Eq. (|T2")) , as in the ID chain the unit cell extends over all p space, while in the 3D condensed matter the cell is 
restricted to x, y € [-R, R]. 

Similar to AEk, in the local approximation, the change in exchange-correlation energy per unit cell is 



A£ cxc (i?) 



\z\<a/2; \x\,\y\<R 



drn 3B (r)e cxc (n w ) - / drn(r) e oxc (n) , 

J \z\<a/2 



(74) 



where e cxc (n) is the exchange-correlation energy (per electron) at density n (see Sec. IHI Al . 

Combining the Coulomb energy, the kinetic energy, and the exchange-correlation energy, the total change in the 
energy per unit cell when 3D condensed matter is formed from ID infinite chains can be written 



where 



AE 3D {R) = AScoui + AE K + AE C 



AE Coul (R) = 8E nD + ^Y1 E QQ^ok) ■ 
(ijk) 



(75) 



(76) 



We calculate AE 3 d (R) as a function of R and locate the minimum to determine the equilibrium chain-chain separation 
2R and the equilibrium energy of the 3D condensed matter. Our method for evaluating various integrals is described 
in the appendix. 
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FIG. 16: The occupation numbers of each m level of infinite Fe chains, for various magnetic field strengths. For B12 > 100, 
only the v = bands are occupied by the electrons (upper panel). For B12 = 5 and 10 the m levels with v = 1 are shown with 
points as well as lines, since there are only a few such occupied levels (lower panel). 
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TABLE V: The energy difference (per unit cell) between the 3D condensed matter and ID chain, AE S = E s — E^, for carbon 
and iron over a range of magnetic field strengths. Energies are given in units of eV for C and keV for Fe. The equilibrium 
chain-chain separation is 27? (in units of the Bohr radius oq). 



B 12 


C 

AE S R 
(eV) 


Fe 

AE S R 
(keV) 


1 


-30 0.200 




5 


-40 0.110 


-0.6 0.150 


10 


-20 0.094 


-0.6 0.115 


100 


-20 0.041 


-2.2 0.054 


500 


-30 0.022 


-2.1 0.025 


1000 


-10 0.017 


-1.3 0.021 



B. Results: 3D condensed matter 

Table |V] presents our numerical results for the equilibrium chain-chain separation 2R = 2R eq and the energy 
difference (per cell) between the 3D condensed matter and ID chain, AE S = E s — Eoo = AEsr>(R = R oq ), for C and 
Fe at various magnetic field strengths. A typical energy curve is shown in Fig. 1171 We see that it is important to 
include the kinetic energy contribution AEk to the 3D energy; without AEk, the energy curve would not have a 
local minimum at a finite R. 

A comparison of the R values in Table |V] with various iron chain electron densities in Fig. [T^] shows that our 
assumption of small electron density overlap between chains is indeed a good approximation. The electron densities 
are slowly-varying at the overlapping region, so using the local (Thomas-Fermi) model to calculate the kinetic energy 
difference is also consistent with the results of our model. Our equilibrium R is within about 15% of the value 
predicted in the uniform cylinder model [see Eq. ((3])] . 

Given our results for AE S and the cohesive energy of ID chains, Qoo = E a — Eoo, we can obtain the cohesive energy 
of 3D condensed matter from 

Q s = E a -E s =E a - (E^ + AE S ) =Q oc - AE S . (77) 

For H and He, we find that |A75 S | is small compared to Qoo and thus Q s ~ Qoo- Figure [TBI depicts Q s and Qoo as a 
function of B for H, He, C, and Fe. 

The only previous quantitative calculation of 3D condensed matter is that by Jones [32j , who finds cohesive energies 
of Q s = 0.60, 0.92 keV for iron at B12 = 5, 10. At these field strengths, our calculation (see Tables HVl and W\ gives 
Q s =E a -E s =Qoo- AE S = 0.08 + 0.6 ~ 0.7 keV and 0.15 + 0.6 ~ 0.75 keV, respectively. 

Note that our calculations and the results presented here assume that the ion spacing along the magnetic axis in 
3D condensed matter, a, is the same as in the ID chain. We have found that if both a and R are allowed to vary, 
the 3D condensed matter energy can be lowered slightly. This correction is most important for relatively low field 
strengths. For example, in the case of Fe at B12 = 10, if we increase a from the ID chain value by 10%, then Qoo 
decreases by about 50 eV, but |A£J S | increases by about 200 eV, so that Q s is increased to ~ 0.9 keV. Given the 
approximate nature of our 3D calculations, we do not explore such refinement in detail in this paper. 

VI. DISCUSSIONS 

Using density functional theory, we have carried out extensive calculations of the cohesive properties of ID infinite 
chains and 3D zero-pressure condensed matter in strong magnetic fields. Our results, presented in various tables, 
figures, and fitting formulae, show that hydrogen, helium, and carbon infinite chains are all bound relative to individual 
atoms for magnetic fields B > 10 12 G, but iron chains are not (significantly) bound until around B ~ 10 14 G. For a 
given zero-pressure condensed matter system, the cohesion along the magnetic axis (chain axis) dominates over chain- 
chain interactions across the magnetic axis at sufficiently strong magnetic fields. But for relative low field strengths 
(e.g. Fe at B < 10 14 G and C at B < a few x 10 12 G), chain-chain interactions play an important role in the cohesion 
of 3D condensed matter. Our calculations show that for the field strengths considered in this paper [B > 10 12 G), 
3D condensed H, He, C and Fe are all bound relative to individual atoms: For C, the cohesive energy Q s = E a — E c 
ranges from - 50 eV at B = 10 12 G to 20 keV at 10 15 G; for Fe, Q s ranges from - 0.8 keV at 10 13 G to 33 keV at 
10 15 G. 
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FIG. 17: The energy (per cell) of 3D condensed matter relative to ID chain as a function of R, for carbon at B12 = 1. The 
nearest-neighbor chain-chain separation in the 3D condensed matter is 2R. The solid curve gives AEso(R) [Eq. (|75|l ] and the 
dashed curve gives only the Coulomb energy AiScoui [Eq. (|76|l ]. 



Our result for the ID infinite chain energy (per cell), E^, is consistent with the energies of finite molecules obtained 
in paper I Q, where we showed that the binding energy (per atom) of the molecule, \Ejf\/N (where En is the ground- 
state energy and N is the number of atoms in the molecule), increases with increasing N, and asymptotes to a constant 
value. The values of \Ejy\/N for various molecules obtained in Ref. Q are always less than |£^oo|. Since the electron 
energy levels in a finite molecule and those in an infinite chain are quite different (the former has discrete states while 
the latter has band structure), and the computations involved are also different, the consistency between the finite 
molecule results and ID chain results provides an important check for the validity of our calculations. 

It is not straightforward to assess the accuracy of our density-functional-theory calculations of infinite chains 
compared to the Hartree-Fock method. For finite molecules with small number of electrons, using the available 
Hartree-Fock results, we have found that density functional theory tends to overestimate the binding energy by 
about 10%, although this does not translate into an appreciable error in the molecular dissociation energy For 
infinite chains, the only previous calculation using the Hartree-Fock method (30j adopted an approximate treatment 
for the electron band structure (e.g., assuming that the electron energy increases as fc 2 /2 as the Bloch wave number 
k increases), which, as we showed in this paper (Sec. IIII C|) . likely resulted in appreciable error to the total chain 
energy. Since the cohesive energy Qoc of the chain involves the difference in the binding energy the ID chain and the 
atom, and because of the statistical nature of density functional theory, we expect that our result for is more 
accurate for heavy elements (C and Fe) than for light elements (H and He). We note that it is very difficult (perhaps 
impractical) to carry out ab initio Hartree-Fock calculations of infinite chains if no approximation is made about the 
electron band structure. This is especially the case in the superstrong magnetic field regime where many Landau 
orbitals are populated. For example, for the Fe chain at B — 10 15 G, one must be dealing with 130 Landau orbitals 
(see Table HV|) . each with its own band structure — this would be a formidable task for any Hartree-Fock calculation. 

We also note that our conclusion about 3D condensed matter is not based on fully self-consistent calculations and 
uses several approximations (Sec. [Vj. Although we have argued that the approximations we adopted are valid and 
our calculation gave reasonable values for the relative binding energies between ID chains and 3D condensed matter, 
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FIG. 18: The cohesive energy as a function of B, for H (dotted line) and He (short-dashed line) infinite chains and C (long-dashed 
lines) and Fe (solid lines) infinite chains (lighter lines) and 3D condensed matter (heavier lines). 

it would be desirable to carry out more definitive calculations of 3D condensed matter. 

Our computed binding energies and equilibrium ion separations of infinite chains and condensed matter agree 
approximately with the simple scaling relations (e.g., Eoo and a as a function of B) derived from the uniform gas 
model (Sec. HIl). We have provided more accurate fitting formulae which will allow one to obtain the cohesive energy at 
various field strengths. Our result for the electron work function (W — |£_f|), however, does not agree with the simple 
scaling relation derived for the uniform electron gas model. For example, we found that W scales more slowly with B 
(7 is significantly smaller than than 2/5) and does not depend strongly on Z (as opposed to the Z A / b dependence for 
the uniform gas model); see Tables ffl llVI This "discrepancy" is understandable since, unlike the B = case, in strong 
magnetic fields the ionization of an atom and binding energy of condensed matter can be very different in values and 
have different dependences on B: for sufficiently large B, the former scales roughly as (In b), while the later scales 
as ~ 6 4 . Our computed electron work function is of order (and usually a fraction of) the ionization energy of the 
corresponding atom, which is generally much smaller than the estimate of W based on uniform gas model. We also 
found that the ionization energy of successively larger (finite) molecules [t| approaches our calculated work function 
for the infinite chain — thus we believe our result for W is reliable. Note that Jones [HJ also found that the work 
function W is almost independent of Z, but his W values scale as B and are much larger than our results for the 
same field strengths. His W values are also larger than the ionization energies of the corresponding atoms. 

Our results for the cohesive energy and work function of condensed matter in strong magnetic fields have significant 
implications for the physical conditions of the outermost layers of magnetized neutron stars and the possible existence 
of "vacuum gap" accelerators in pulsars. We plan to investigate these issues in the future. 
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APPENDIX A: TECHNICAL DETAILS AND NUMERICAL METHOD 



1. Evaluating the integrals in the Kohn-Sham equations 

The most computation-intensive term in the modified Kohn-Sham equations [Eqs. (|35j) and (j4Tj) ] is the direct 
electron-electron interaction term 



^ee,m (^) 



dr i dt 



|z'|<a(jV Q +l/2) 



,\ Wm\ 2 (p)n(r<) 
|r-r'| 



(Al) 



The evaluation of this term is the rate-limiting step in the entire energy calculation. The integral is over four variables 
(p, p' , z', and </> or 6 — <//), so it requires some simplification to become tractable. To simplify the integral we use the 
identity (e.g., Ref. g|) 



-| °° pOO 

—ii= E / dqe m ^-^J n ( q p)J n (qp')e^ z - z ' 



where J n {z) is the nth order Bessel function of the first kind. Then 

n(r') 



V ee {r) 



dr' 



|z'|<a(jV Q + l/2) 
a(N Q + l/2) 



r — r 



/m«Q + l/^ rOO rOO 

dz' dqJ (qp) / p' dp' n{p',z')J {qp') 
-a(N Q + l/2) JO UO 



exp(— q\z — z'\ 



and 



V, ,,„(r) : : / dr ± \W m \ 2 (p)V ee (r) 

r a(N Q + l/2) 
i(JV Q + l/2) 



ra(N Q + l/2) poo r r 

4tt 2 / dz' / / 

J -a(N +l/2) JO Uo 



pdp\W m \ 2 (p)Jo(qp) 



p' dp' n(p',z')J a (qp') 



exp(— q\z — z'\ 



Using Eq. (fT6|) for the electron density distribution, Eq. (|A6| becomes 

ra(N Q + l/2) 



/a(i\ Q + L/. 
dz' f^, v ,{z') I dqG m (q)G m ,(q)exp(-q\z - z'\) . 
,„ „ -o(jv q +i/2) Jo 



where 



G m (q) = 2ir / pdp\W m \ 2 {p) J (qp) 
Jo 

= exp(-< Z 2 /2) J L m (q 2 /2), 



and 



L m (x) = 



e x d n 



■{x m e~ x ) 



to! dx r 

is the Laguerre polynomial of order m. These polynomials can be calculated using the recurrence relation 

mL m (x) = (2m - 1 - x)L m _i(x) - (to - l)L m _ 2 (x) , 



(A2) 

(A3) 
(A4) 

(A5) 

'!)• 

(A6) 

(A7) 

(A8) 
(A9) 

(A10) 

(All) 



with Lq(x) = 1 and L±(x) = 1 — x. 

Using the method outlined above the original four-dimensional integral in Eq. (|A1[) reduces to a two-dimensional 
integral. Once a value for z is specified, the integral can be evaluated using a quadrature algorithm (such as the 
Romberg integration method [48j]). 
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2. Evaluating the integrals in the calculation of 3D condensed matter 



For the 3D condensed matter calculation, we simplify the energy integrals of the nearest-neighbor interactions in a 
way similar to that for the infinite chain calculation. To do this, we require Eq. (|A2|) and one additional identity of 
Bessel functions: 



Jo(qVa 2 + b 2 -2abcos8) = ^ e inS ' J n (qa)J n (qb) . 



n— — oo 



With these equations the ion-electron nearest-neighbor energy term [Eq. ([65]) ] becomes 

n(r) 



-Ze l 



dr ■ 

\<a/2 |r — r nn 

• a/2 



/a / z />oo />oo 

dz / dqJo(2Rq) / pdpn(p, z)Jq 
-a/2 JO Uo 



(qp) 



'-a/ 
r a/2 



exp(— q\z — a/2|) 



Ze 2 J2 / dzf 2 v (z) \ 



dz fmA z ) / dqJ a (2Rq)G m {q)exp{-q\z-a/2\). 



The electron-electron energy term [Eq. (|66|) ] becomes 

1 <" r , n(r)n(r') 



Edh,nn[n] 



dvdr 

\z\<a/2,\z'\<a/2 |r-(r' + r nn )| 



e 2 o 
— 2tt 

2 



a/2 
-a/2 

27T 



dz 
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dz 1 / p' dp' n(p', z') x 



/>Z7T /■ OO 

/ #' / J (g|r' ± +r ±!nn |) 
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-q\z-z' —a/2 



where is the angle of r' ± + r^ nn in the (p, 0, 2;) cylindrical coordinate system 



e 2 /-Q/2 /-o/2 

-Bdir.nnW = — -47T 2 / / . 
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pdpn(p, z)J (qp) 



a/2 



p'dp'n(p',z')J (qp') 
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(A12) 

(A13) 
(A14) 
(A15) 

(A16) 



(A17) 



(A18) 



dzf 2 „(z) dz'f^ v ,{z') dqM2Rq)G m (q)G' m (q)exp(-q\z-z'-a/2\). 



(A19) 



Notice that the infinite chain expression for the nearest-neighbor electron-electron interaction energy is recovered 
when R — and a/2 is replaced by ±a. 



3. Solving the differential equations and the total energy self-consistently 

The Kohn-Sham equations [Eqs. ([3"5]) and (|4"Tj) ] are solved on a grid in z. Because of symmetry we only need to 
consider z > 0, with z — coincident with an ion. The number and spacing of the z grid points determine how 
accurately the equations can be solved. In this paper we have attempted to calculate ground-state chain energies 
to better than 0.1% numerical accuracy. This requires approximately (depending on Z and B) 33 grid points for 
each unit cell and 3 cells (for Nq = 1 there are three cells that require exact treatment: the cell under consideration 
z e [—a/2, a/2] and its nearest neighbors; the rest of the cells enter the calculation only through their quadrupolc 
moments). The grid spacing is chosen to be constant from the center out to the edge of the cell. The shape of the 
wave function is found within one cell and then copied to the other cells. 



31 



For integration with respect to p, p' , or q (e.g., when calculating the direct electron-electron interaction term), our 
goal of 0.1% accuracy for the total energy requires an accuracy of approximately 10 -5 in the integral. A variable-step- 
size integration routine is used for each such integral, where the number of points in the integration grid is increased 
until the error in the integration is within the desired accuracy. 

We discussed the boundary conditions for the wave function solutions to the Kohn-Sham equations (see Sec. IIII Cj) . 
The only other requirement we have for these wave functions is that the magnitude of each wave function has the 
correct number (v) of nodes per cell (see Fig. [2]). In practice, to find fm»o( z ) an d f m v- K / a {z) we integrate Eqs. (f3l)|) 
and JIT]) from one edge of the z grid (e.g., z = a/2) and shoot toward the center (z = 0), adjusting e mv (k = 0) and 
£ m „(7r/o) until the correct boundary condition is satisfied. For the other k values with energies between these two 
extremes, we use the given energy to find a wave function and calculate the k that solves the boundary condition 
Eq. (|5Tj) , as discussed in Sec. IIII Cl 

There are two parts to our procedure for finding / m „fc(z), s mu (k), and <j mv self-consistently: (i) determining 
the longitudinal wave functions f m vk and periodic potential self-consistently, and (ii) determining the electron level 
occupations a mv self-consistently. 

To determine the f m vk wave functions self-consistently, a trial set of wave functions and <7 m „ values is first used to 
calculate the potential as a function of z, and that potential is used to calculate a new set of wave functions. These 
new wave functions are then used to find a new potential, and the process is repeated until consistency is reached. 
In practice, we find that f m vk{z) = works well as the trial wave function and a linear spread of a mv from oq v = 1 
to a nmV = works well for the trial a values. Convergence can be achieved in four or five iterations. To prevent 
overcorrection from one iteration to the next, the actual potential used for each iteration is a combination of the 
newly-generated potential and the old potential from the previous iteration (the weighting used is roughly 30% old, 
70% new). 

To determine the <j mv level occupations self-consistently, we first find the wave functions and eigenvalues e mi/ (k) 
as a function of k self-consistently as described above. With this information, and given a Fermi level energy sf, we 
can calculate new a values, using the equations in Sec. III. The Fermi level energy is adjusted until a mu = Z using 
Newton's method. These new a mv values are used to re-calculate the wave functions self-consistently. This process is 
repeated until self-consistency is reached, which is typically after about three (for hydrogen at 10 12 G) to twelve (for 
iron at 2 x 10 15 G) full iterations. 
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